TODO

  • add distribution plots for all (sub)scales
  • not all subscales in TRT plots
  • add internal consistency subscales
  • add bar plot for IC
  • add second set of correlations for t2? could put in supplementaries
  • could include the PTTA and PTT analyses very briefly as exploratory, ie just to replicate the weird claim that spatial and emotional PT are correlated (no) and that the PTT measures PT (well, doesn’t correlate with other PT tasks).

Overview

All included measures are used as measures of empathy or to make conclusions about empathy. While some may be used to examine different facets of empathy, this distinction is a) not theoretically agreed upon, as there is little consensus around the structure of empathy, and therefore b) not necessarily reflected in their measurement (i.e., low correlations may represent poor convergent validity rather than good discriminant validity).

This study examines a) the reliability of these scales (internal consistency, test-retest reliability) relative to the original published estimates, and b) the risk of jingle fallacy (via correlations among [sub]scales).

NB EQ subscales based on Lawrence et al., 2004.

Dependencies

library(readr)
library(dplyr)
library(tidyr)
library(stringr)
library(forcats)
library(irr)
library(psych)
library(lavaan)
library(janitor)
library(ggplot2)
library(patchwork)
library(ggstance)
library(ggrain)
library(scales)
library(mirt)
library(correlation)
library(see)
library(insight)
library(purrr)
library(knitr)
library(kableExtra)

options(scipen=999)

dir.create("plots")

Load data

data_processed <- read_csv("../data/processed/data_processed.csv") |>
  select(-starts_with("PTP_") & -starts_with("PTTA_"))

data_processed_after_exclusions <- data_processed |>
  filter(exclude_master == "include")

Demographics

Sample sizes and exclusions

data_processed |>
  count(timepoint, exclude_master) |>
  kable(align = "r") |>
  kable_classic(full_width = FALSE)
timepoint exclude_master n
1 exclude 2
1 include 134
2 exclude 2
2 include 99

Analytic sample (timepoint 1)

data_processed_after_exclusions |>
  filter(timepoint == 1) |>
  count(gender) |>
  kable(align = "r") |>
  kable_classic(full_width = FALSE)
gender n
Non binary 1
female 91
male 42 
data_processed_after_exclusions |>
  filter(timepoint == 1) |>
  dplyr::summarize(mean_age = mean(age),
                   sd_age = sd(age),
                   min_age = min(age),
                   max_age = max(age)) |>
  pivot_longer(cols = everything()) |>
  mutate(across(where(is.numeric), ~ round_half_up(.x, digits = 2))) |>
  kable(align = "r") |>
  kable_classic(full_width = FALSE)
name value
mean_age 38.23
sd_age 13.04
min_age 19.00
max_age 75.00 
data_processed_after_exclusions |>
  distinct(unique_id, .keep_all = TRUE) |>
  ggplot(aes(age)) +
  geom_histogram(binwidth = 5) +
  theme_linedraw()

Analyses

Test-retest

\(ICC_{2,1}\)

Results

data_reshaped <- data_processed_after_exclusions |>
  select(unique_id, timepoint, SITES, IRI_FS, IRI_EC, IRI_PT, IRI_PD, IRI_total, EQ_total, EQ_cog, EQ_aff, EQ_soc, PET_total, EYES_total) |>
  pivot_wider(id_cols = unique_id,
              names_from = timepoint,
              values_from = c(SITES, IRI_FS, IRI_EC, IRI_PT, IRI_PD, IRI_total, EQ_total, EQ_cog, EQ_aff, EQ_soc, PET_total, EYES_total))

tidy_icc_2_1 <- function(data, scale){
  res <- data |>
    select(starts_with(scale)) |>
    irr::icc(model = "twoway", type = "agreement", unit = "single",
             r0 = 0, conf.level = 0.95)
  
  tibble(scale = scale,
         icc_2.1 = res$value,
         ci_lower = res$lbound,
         ci_upper = res$ubound)
}

res <- bind_rows(
  tidy_icc_2_1(data_reshaped, "SITES"),
  tidy_icc_2_1(data_reshaped, "IRI_FS"),
  tidy_icc_2_1(data_reshaped, "IRI_EC"),
  tidy_icc_2_1(data_reshaped, "IRI_PT"),
  tidy_icc_2_1(data_reshaped, "IRI_PD"),
  tidy_icc_2_1(data_reshaped, "IRI_total"),
  tidy_icc_2_1(data_reshaped, "EQ_total"),
  tidy_icc_2_1(data_reshaped, "EQ_cog"),
  tidy_icc_2_1(data_reshaped, "EQ_aff"),
  tidy_icc_2_1(data_reshaped, "EQ_soc"),
  tidy_icc_2_1(data_reshaped, "PET_total"),
  tidy_icc_2_1(data_reshaped, "EYES_total")
)
  
res |>
  mutate(across(where(is.numeric), ~ round_half_up(.x, digits = 2))) |>
  kable(align = "r") |>
  kable_classic(full_width = FALSE)
scale icc_2.1 ci_lower ci_upper
SITES 0.77 0.67 0.84
IRI_FS 0.88 0.82 0.92
IRI_EC 0.87 0.82 0.91
IRI_PT 0.81 0.73 0.87
IRI_PD 0.88 0.82 0.92
IRI_total 0.90 0.86 0.93
EQ_total 0.86 0.80 0.91
EQ_cog 0.84 0.76 0.89
EQ_aff 0.77 0.68 0.84
EQ_soc 0.74 0.64 0.82
PET_total 0.86 0.80 0.91
EYES_total 0.73 0.63 0.81

Plots

Bar

p_trt_scales <- res |>
  filter(str_detect(scale, "SITES") | str_detect(scale, "_total")) |>
  mutate(scale = str_remove(scale, "_total")) |>
  ggplot(aes(icc_2.1, fct_rev(scale))) +
  geom_vline(xintercept = 0, linetype = "dotted") +
  geom_linerangeh(aes(xmin = ci_lower, xmax = ci_upper), position = position_dodge(width = .3)) +
  geom_point(position = position_dodge(width = .3)) +
  scale_x_continuous(breaks = breaks_pretty(n = 10),
                     limits = c(-.2, 1),
                     name = bquote(ICC["2,1"])) +
  theme_linedraw() +
  ylab("Scale")

p_trt_subscales <- res |>
  filter(!str_detect(scale, "_total") & !str_detect(scale, "SITES")) |>
  mutate(scale = str_replace(scale, "_", " ")) |>
  mutate(scale = str_replace(scale, "doublereversed", "double")) |>
  ggplot(aes(icc_2.1, fct_rev(scale))) +
  geom_vline(xintercept = 0, linetype = "dotted") +
  geom_linerangeh(aes(xmin = ci_lower, xmax = ci_upper), position = position_dodge(width = .3)) +
  geom_point(position = position_dodge(width = .3)) +
  scale_x_continuous(breaks = breaks_pretty(n = 10),
                     limits = c(-.2, 1),
                     name = bquote(ICC["2,1"])) +
  theme_linedraw() +
  ylab("Subscale")

p_trt_scales + 
  p_trt_subscales + 
  plot_layout(ncol = 1, heights = c(1, 1.35))

ggsave("plots/trt.png", 
       width = 6,
       height = 5,
       dpi = 600)

Raincloud

round_half_up_trailing <- function(x, digits = 2) {
  x_rounded <- round_half_up(x, digits = digits)
  sprintf("%.*f", digits, x_rounded)
}

res_print <- res |>
  mutate(across(where(is.numeric), ~ round_half_up_trailing(.x, digits = 2))) |>
  mutate(string = paste0("= ", icc_2.1, ", 95% CI [", ci_lower, ", ", ci_upper, "]")) |>
  select(scale, string)

data_reshaped_plotting <- data_processed_after_exclusions |>
  select(unique_id, timepoint, SITES, IRI_FS, IRI_EC, IRI_PT, IRI_PD, IRI_total, EQ_total, EQ_cog, EQ_aff, EQ_soc, PET_total, EYES_total) |>
  pivot_longer(cols = c(SITES, IRI_FS, IRI_EC, IRI_PT, IRI_PD, IRI_total, EQ_total, EQ_cog, EQ_aff, EQ_soc, PET_total, EYES_total),
               names_to = "scale",
               values_to = "score")

p_trt_1 <- data_reshaped_plotting |>
  filter(scale == "IRI_total") |>
  ggplot(aes(as.factor(timepoint), score, fill = as.factor(timepoint))) +
  geom_rain(alpha = .5, rain.side = 'f1x1', id.long.var = "unique_id") +
  theme_classic() +
  scale_fill_manual(values=c("dodgerblue", "darkorange")) +
  guides(fill = 'none', color = 'none') +
  ylab("Score") +
  xlab("Time point") +
  ggtitle(bquote("Interpersonal Reactivity Index:" ~ ICC["2,1"] ~ .(res_print |> filter(scale == "IRI_total") |> pull(string)))) +
  scale_y_continuous(breaks = scales::breaks_pretty(n = 6))

p_trt_2 <- data_reshaped_plotting |>
  filter(scale == "IRI_PT") |>
  ggplot(aes(as.factor(timepoint), score, fill = as.factor(timepoint))) +
  geom_rain(alpha = .5, rain.side = 'f1x1', id.long.var = "unique_id") +
  theme_classic() +
  scale_fill_manual(values=c("dodgerblue", "darkorange")) +
  guides(fill = 'none', color = 'none') +
  ylab("Score") +
  xlab("Time point") +
  ggtitle(bquote("IRI - Perspective Taking Subscale:" ~ ICC["2,1"] ~ .(res_print |> filter(scale == "IRI_PT") |> pull(string)))) +
  scale_y_continuous(breaks = scales::breaks_pretty(n = 6))

p_trt_3 <- data_reshaped_plotting |>
  filter(scale == "SITES") |>
  ggplot(aes(as.factor(timepoint), score, fill = as.factor(timepoint))) +
  geom_rain(alpha = .5, rain.side = 'f1x1', id.long.var = "unique_id", likert = TRUE) +
  theme_classic() +
  scale_fill_manual(values=c("dodgerblue", "darkorange")) +
  guides(fill = 'none', color = 'none') +
  ylab("Score") +
  xlab("Time point") +
  ggtitle(bquote("Single Item Trait Empathy Scale:" ~ ICC["2,1"] ~ .(res_print |> filter(scale == "SITES") |> pull(string)))) +
  scale_y_continuous(breaks = scales::breaks_pretty(n = 6))

p_trt_4 <- data_reshaped_plotting |>
  filter(scale == "PET_total") |>
  ggplot(aes(as.factor(timepoint), score, fill = as.factor(timepoint))) +
  geom_rain(alpha = .5, rain.side = 'f1x1', id.long.var = "unique_id") +
  theme_classic() +
  scale_fill_manual(values=c("dodgerblue", "darkorange")) +
  guides(fill = 'none', color = 'none') +
  ylab("Score") +
  xlab("Time point") +
  ggtitle(bquote("Pictoral Empathy Test:" ~ ICC["2,1"] ~ .(res_print |> filter(scale == "PET_total") |> pull(string)))) +
  scale_y_continuous(breaks = scales::breaks_pretty(n = 6))

p_trt_5 <- data_reshaped_plotting |>
  filter(scale == "EQ_total") |>
  ggplot(aes(as.factor(timepoint), score, fill = as.factor(timepoint))) +
  geom_rain(alpha = .5, rain.side = 'f1x1', id.long.var = "unique_id") +
  theme_classic() +
  scale_fill_manual(values=c("dodgerblue", "darkorange")) +
  guides(fill = 'none', color = 'none') +
  ylab("Score") +
  xlab("Time point") +
  ggtitle(bquote("Emotional Quotient:" ~ ICC["2,1"] ~ .(res_print |> filter(scale == "EQ_total") |> pull(string)))) +
  scale_y_continuous(breaks = scales::breaks_pretty(n = 6))

p_trt_6 <- data_reshaped_plotting |>
  filter(scale == "EYES_total") |>
  ggplot(aes(as.factor(timepoint), score, fill = as.factor(timepoint))) +
  geom_rain(alpha = .5, rain.side = 'f1x1', id.long.var = "unique_id") +
  theme_classic() +
  scale_fill_manual(values=c("dodgerblue", "darkorange")) +
  guides(fill = 'none', color = 'none') +
  ylab("Score") +
  xlab("Time point") +
  ggtitle(bquote("Reading the Mind in the Eyes Test:" ~ ICC["2,1"] ~ .(res_print |> filter(scale == "EYES_total") |> pull(string)))) +
  scale_y_continuous(breaks = scales::breaks_pretty(n = 6))

p_trt_1 + p_trt_2 + p_trt_3 + p_trt_4 +
  p_trt_5 + p_trt_6 + 
  plot_layout(ncol = 2)

Internal consistency

Cronbach’s \(\alpha\)

data_subset <- data_processed_after_exclusions |>
  dplyr::select(timepoint, 
                starts_with("IRI_"), 
                starts_with("EQ_"),
                starts_with("PET_"),
                starts_with("EYES_")) |>
  dplyr::select(-IRI_FS, -IRI_EC, -IRI_PT, -IRI_PD, -EQ_cog, -EQ_aff, -EQ_soc) |>
  dplyr::select(!ends_with("_completeness") & !ends_with("_total"))

tidy_alpha <- function(data, time, scale){
  res <- data |>
    filter(timepoint == time) |>
    dplyr::select(starts_with(scale)) |>
    psych::alpha()
  
  res2 <- res$feldt 
  
  tibble(scale = scale,
         timepoint = time,
         alpha = res2$alpha$raw_alpha,
         ci_lower = res2$lower.ci$raw_alpha,
         ci_upper = res2$upper.ci$raw_alpha)
}

res <- bind_rows(
  tidy_alpha(data = data_subset, time = 1, scale = "IRI_"),
  #tidy_alpha(data = data_subset, time = 2, scale = "IRI_"),
  tidy_alpha(data = data_subset, time = 1, scale = "EQ_"),
  #tidy_alpha(data = data_subset, time = 2, scale = "EQ_"),
  tidy_alpha(data = data_subset, time = 1, scale = "PET_"),
  #tidy_alpha(data = data_subset, time = 2, scale = "PET_"),
  tidy_alpha(data = data_subset, time = 1, scale = "EYES_"),
  #tidy_alpha(data = data_subset, time = 2, scale = "EYES_")
)
## Some items ( EQ_57 ) were negatively correlated with the first principal component and 
## probably should be reversed.  
## To do this, run the function again with the 'check.keys=TRUE' option
res |>
  mutate(across(where(is.numeric), ~ round_half_up(.x, digits = 2))) |>
  kable(align = "r") |>
  kable_classic(full_width = FALSE)
scale timepoint alpha ci_lower ci_upper
IRI_ 1 0.88 0.85 0.91
EQ_ 1 0.89 0.87 0.92
PET_ 1 0.92 0.89 0.94
EYES_ 1 0.68 0.60 0.76

Correlations

T1

dat <- data_processed |>
  filter(timepoint == 1)

Scale correlations

scale_cors <- dat |>
  select(SITES, ends_with("_total")) |>
  rename_with(~ str_remove(.x, "_total")) %>%
  select(order(colnames(.))) |>
  correlation()

res <- scale_cors |>
  as.data.frame() |>
  select(Parameter1, Parameter2, r, ci_low = CI_low, ci_high = CI_high, p) 

res |>
  mutate(p = insight::format_p(p)) |>
  mutate(across(where(is.numeric), ~ round_half_up(.x, digits = 2))) |>
  kable(align = "r") |>
  kable_classic(full_width = FALSE)
Parameter1 Parameter2 r ci_low ci_high p
EQ EYES 0.26 0.09 0.41 p = 0.011
EQ IRI 0.49 0.36 0.61 p < .001
EQ PET 0.42 0.27 0.55 p < .001
EQ SITES 0.54 0.41 0.65 p < .001
EYES IRI 0.09 -0.08 0.25 p = 0.635
EYES PET 0.11 -0.06 0.27 p = 0.628
EYES SITES 0.07 -0.10 0.24 p = 0.635
IRI PET 0.65 0.55 0.74 p < .001
IRI SITES 0.57 0.44 0.67 p < .001
PET SITES 0.44 0.30 0.57 p < .001 
p1_t1 <- 
  ggplot(res, aes(r)) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  geom_histogram(binwidth = 0.1, boundary = 0) +
  scale_x_continuous(limits = c(-1, 1), 
                     name = "Correlation", 
                     breaks = scales::breaks_pretty(n = 10)) + 
  scale_y_continuous(name = "Frequency", 
                     breaks = scales::breaks_pretty(n = 7)) + 
  theme_linedraw() +
  theme(panel.grid.minor.y = element_blank())

#p1
  
p2_t1 <- 
  scale_cors |>
  summary(redundant = TRUE) |>
  plot(show_data = "points",
       stars = FALSE) +
  theme_linedraw() +
  ggtitle("")

#p2 + coord_fixed(ratio = 1)

# scale_cors |>
#   summary() |>
#   plot(show_data = "points",
#        stars = FALSE) +
#   theme_linedraw() +
#   ggtitle("") +
#   coord_fixed(ratio = 1)

average_cors_all <- res |>
  dplyr::summarize(mean_r = mean(r),
                   sd_r = sd(r),
                   median_r = median(r),
                   mad_r = mad(r),
                   min_r = min(r),
                   max_r = max(r)) 

average_cors_all |>
  round_half_up(digits = 2) |>
  kable() |>
  kable_classic(full_width = FALSE)
mean_r sd_r median_r mad_r min_r max_r
0.36 0.22 0.43 0.23 0.07 0.65 
p1_t1 + ggtitle("Correlation histogram") +
  p2_t1 + ggtitle("Correlation matrix") + 
  plot_layout(ncol = 1, heights = c(.4, 1))

ggsave("plots/cor_scales_t1.png", 
       width = 6,
       height = 7,
       dpi = 600)

Subscale correlations

subscale_cors <- dat |>
  select(SITES, 
         IRI_FS, IRI_EC, IRI_PT, IRI_PD, 
         EQ_cog, EQ_aff, EQ_soc, 
         PET_total,
         EYES_total) |>
  rename_with(~ str_remove(.x, "_total")) |>
  correlation() 

subscale_cors |>
  as.data.frame() |>
  select(Parameter1, Parameter2, r, ci_low = CI_low, ci_high = CI_high, p) |>
  mutate(p = insight::format_p(p)) |>
  mutate(across(where(is.numeric), ~ round_half_up(.x, digits = 2))) |>
  kable(align = "r") |>
  kable_classic(full_width = FALSE)
Parameter1 Parameter2 r ci_low ci_high p
SITES IRI_FS 0.33 0.17 0.48 p = 0.002
SITES IRI_EC 0.68 0.58 0.76 p < .001
SITES IRI_PT 0.54 0.41 0.65 p < .001
SITES IRI_PD 0.12 -0.05 0.28 p > .999
SITES EQ_cog 0.48 0.34 0.60 p < .001
SITES EQ_aff 0.55 0.42 0.66 p < .001
SITES EQ_soc 0.29 0.13 0.44 p = 0.014
SITES PET 0.44 0.30 0.57 p < .001
SITES EYES 0.07 -0.10 0.24 p > .999
IRI_FS IRI_EC 0.52 0.38 0.63 p < .001
IRI_FS IRI_PT 0.36 0.20 0.50 p < .001
IRI_FS IRI_PD 0.28 0.12 0.43 p = 0.018
IRI_FS EQ_cog 0.16 -0.01 0.32 p = 0.771
IRI_FS EQ_aff 0.47 0.33 0.59 p < .001
IRI_FS EQ_soc 0.16 -0.01 0.32 p = 0.771
IRI_FS PET 0.50 0.36 0.62 p < .001
IRI_FS EYES 0.10 -0.07 0.27 p > .999
IRI_EC IRI_PT 0.65 0.54 0.74 p < .001
IRI_EC IRI_PD 0.25 0.09 0.40 p = 0.057
IRI_EC EQ_cog 0.44 0.30 0.57 p < .001
IRI_EC EQ_aff 0.66 0.55 0.74 p < .001
IRI_EC EQ_soc 0.26 0.09 0.41 p = 0.045
IRI_EC PET 0.61 0.49 0.71 p < .001
IRI_EC EYES 0.05 -0.12 0.22 p > .999
IRI_PT IRI_PD -0.01 -0.18 0.16 p > .999
IRI_PT EQ_cog 0.44 0.30 0.57 p < .001
IRI_PT EQ_aff 0.55 0.42 0.66 p < .001
IRI_PT EQ_soc 0.31 0.14 0.45 p = 0.007
IRI_PT PET 0.43 0.28 0.56 p < .001
IRI_PT EYES 0.08 -0.09 0.25 p > .999
IRI_PD EQ_cog -0.18 -0.34 -0.01 p = 0.488
IRI_PD EQ_aff 0.21 0.04 0.37 p = 0.191
IRI_PD EQ_soc 0.17 0.00 0.33 p = 0.569
IRI_PD PET 0.32 0.16 0.47 p = 0.003
IRI_PD EYES 0.01 -0.16 0.18 p > .999
EQ_cog EQ_aff 0.51 0.37 0.63 p < .001
EQ_cog EQ_soc 0.27 0.10 0.42 p = 0.034
EQ_cog PET 0.25 0.08 0.40 p = 0.061
EQ_cog EYES 0.15 -0.02 0.31 p = 0.771
EQ_aff EQ_soc 0.53 0.40 0.64 p < .001
EQ_aff PET 0.56 0.43 0.66 p < .001
EQ_aff EYES 0.25 0.08 0.40 p = 0.061
EQ_soc PET 0.22 0.05 0.37 p = 0.156
EQ_soc EYES 0.33 0.17 0.47 p = 0.003
PET EYES 0.11 -0.06 0.27 p > .999 
p3_t1 <- 
  ggplot(subscale_cors, aes(r)) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  geom_histogram(binwidth = 0.1, boundary = 0) +
  scale_x_continuous(limits = c(-1, 1), 
                     name = "Correlation", 
                     breaks = scales::breaks_pretty(n = 10)) + 
  scale_y_continuous(name = "Frequency", 
                     breaks = scales::breaks_pretty(n = 7)) + 
  theme_linedraw() +
  theme(panel.grid.minor.y = element_blank())

#p3_t1

# code repeated in order to use the rename, which screws up the table but is needed for the plot
p4_t1 <- dat |>
  select(SITES, 
         IRI_FS, IRI_EC, IRI_PT, IRI_PD, 
         EQ_cog, EQ_aff, EQ_soc, 
         PET_total,
         EYES_total) |>
  rename_with(~ str_remove(.x, "_total")) |>
  rename_with(~ str_replace(.x, "_", "\n")) %>%
  select(order(colnames(.))) |>
  correlation() %>%
  mutate(across(where(is.numeric), ~ round_half_up(., 2))) |>
  summary(redundant = TRUE) |>
  plot(show_data = "points",
       stars = FALSE) +
  theme_linedraw() +
  ggtitle("")

#p4_t1
p3_t1 + ggtitle("Correlation histogram") +
  p4_t1 + ggtitle("Correlation matrix") + 
  plot_layout(ncol = 1, heights = c(.4, 1))

ggsave("plots/cor_subscales_t1.png", 
       width = 7,
       height = 9,
       dpi = 600)

T2

dat <- data_processed |>
  filter(timepoint == 2)

Scale correlations

scale_cors <- dat |>
  select(SITES, ends_with("_total")) |>
  rename_with(~ str_remove(.x, "_total")) %>%
  select(order(colnames(.))) |>
  correlation()

res <- scale_cors |>
  as.data.frame() |>
  select(Parameter1, Parameter2, r, ci_low = CI_low, ci_high = CI_high, p) 

res |>
  mutate(p = insight::format_p(p)) |>
  mutate(across(where(is.numeric), ~ round_half_up(.x, digits = 2))) |>
  kable(align = "r") |>
  kable_classic(full_width = FALSE)
Parameter1 Parameter2 r ci_low ci_high p
EQ EYES 0.27 0.08 0.44 p = 0.028
EQ IRI 0.47 0.30 0.61 p < .001
EQ PET 0.32 0.13 0.49 p = 0.007
EQ SITES 0.52 0.36 0.65 p < .001
EYES IRI 0.12 -0.08 0.30 p = 0.665
EYES PET 0.11 -0.09 0.30 p = 0.665
EYES SITES 0.12 -0.07 0.31 p = 0.665
IRI PET 0.59 0.45 0.71 p < .001
IRI SITES 0.50 0.34 0.63 p < .001
PET SITES 0.27 0.08 0.45 p = 0.028 
p1_t2 <- 
  ggplot(res, aes(r)) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  geom_histogram(binwidth = 0.1, boundary = 0) +
  scale_x_continuous(limits = c(-1, 1), 
                     name = "Correlation", 
                     breaks = scales::breaks_pretty(n = 10)) + 
  scale_y_continuous(name = "Frequency", 
                     breaks = scales::breaks_pretty(n = 7)) + 
  theme_linedraw() +
  theme(panel.grid.minor.y = element_blank())

#p1_t2
  
p2_t2 <- 
  scale_cors |>
  summary(redundant = TRUE) |>
  plot(show_data = "points",
       stars = FALSE) +
  theme_linedraw() +
  ggtitle("")

#p2_t2 + coord_fixed(ratio = 1)

# scale_cors |>
#   summary() |>
#   plot(show_data = "points",
#        stars = FALSE) +
#   theme_linedraw() +
#   ggtitle("") +
#   coord_fixed(ratio = 1)

average_cors_all <- res |>
  dplyr::summarize(mean_r = mean(r),
                   sd_r = sd(r),
                   median_r = median(r),
                   mad_r = mad(r),
                   min_r = min(r),
                   max_r = max(r)) 

average_cors_all |>
  round_half_up(digits = 2) |>
  kable() |>
  kable_classic(full_width = FALSE)
mean_r sd_r median_r mad_r min_r max_r
0.33 0.18 0.3 0.26 0.11 0.59 
p1_t2 + ggtitle("Correlation histogram") +
  p2_t2 + ggtitle("Correlation matrix") + 
  plot_layout(ncol = 1, heights = c(.4, 1))

ggsave("plots/cor_scales_t2.png", 
       width = 6,
       height = 7,
       dpi = 600)

Subscale correlations

subscale_cors <- dat |>
  select(SITES, 
         IRI_FS, IRI_EC, IRI_PT, IRI_PD, 
         EQ_cog, EQ_aff, EQ_soc, 
         PET_total,
         EYES_total) |>
  rename_with(~ str_remove(.x, "_total")) |>
  correlation() 

subscale_cors |>
  as.data.frame() |>
  select(Parameter1, Parameter2, r, ci_low = CI_low, ci_high = CI_high, p) |>
  mutate(p = insight::format_p(p)) |>
  mutate(across(where(is.numeric), ~ round_half_up(.x, digits = 2))) |>
  kable(align = "r") |>
  kable_classic(full_width = FALSE)
Parameter1 Parameter2 r ci_low ci_high p
SITES IRI_FS 0.20 0.00 0.38 p = 0.970
SITES IRI_EC 0.65 0.51 0.75 p < .001
SITES IRI_PT 0.48 0.31 0.62 p < .001
SITES IRI_PD 0.09 -0.11 0.28 p > .999
SITES EQ_cog 0.45 0.27 0.59 p < .001
SITES EQ_aff 0.48 0.31 0.62 p < .001
SITES EQ_soc 0.13 -0.07 0.32 p > .999
SITES PET 0.27 0.08 0.45 p = 0.146
SITES EYES 0.12 -0.07 0.31 p > .999
IRI_FS IRI_EC 0.49 0.32 0.62 p < .001
IRI_FS IRI_PT 0.29 0.10 0.46 p = 0.086
IRI_FS IRI_PD 0.20 0.00 0.38 p = 0.970
IRI_FS EQ_cog 0.12 -0.08 0.31 p > .999
IRI_FS EQ_aff 0.47 0.30 0.61 p < .001
IRI_FS EQ_soc 0.01 -0.19 0.20 p > .999
IRI_FS PET 0.46 0.29 0.60 p < .001
IRI_FS EYES 0.13 -0.07 0.31 p > .999
IRI_EC IRI_PT 0.56 0.40 0.68 p < .001
IRI_EC IRI_PD 0.12 -0.08 0.31 p > .999
IRI_EC EQ_cog 0.39 0.21 0.55 p = 0.002
IRI_EC EQ_aff 0.71 0.59 0.79 p < .001
IRI_EC EQ_soc 0.10 -0.10 0.29 p > .999
IRI_EC PET 0.55 0.40 0.68 p < .001
IRI_EC EYES 0.02 -0.18 0.21 p > .999
IRI_PT IRI_PD -0.05 -0.24 0.15 p > .999
IRI_PT EQ_cog 0.50 0.34 0.64 p < .001
IRI_PT EQ_aff 0.51 0.35 0.64 p < .001
IRI_PT EQ_soc 0.14 -0.06 0.33 p > .999
IRI_PT PET 0.32 0.13 0.49 p = 0.035
IRI_PT EYES 0.12 -0.08 0.31 p > .999
IRI_PD EQ_cog -0.22 -0.40 -0.02 p = 0.643
IRI_PD EQ_aff 0.10 -0.10 0.29 p > .999
IRI_PD EQ_soc 0.00 -0.20 0.19 p > .999
IRI_PD PET 0.25 0.06 0.43 p = 0.286
IRI_PD EYES 0.05 -0.15 0.24 p > .999
EQ_cog EQ_aff 0.48 0.31 0.62 p < .001
EQ_cog EQ_soc 0.23 0.03 0.40 p = 0.555
EQ_cog PET 0.11 -0.09 0.30 p > .999
EQ_cog EYES 0.20 0.00 0.38 p = 0.970
EQ_aff EQ_soc 0.34 0.15 0.50 p = 0.018
EQ_aff PET 0.47 0.30 0.61 p < .001
EQ_aff EYES 0.27 0.08 0.45 p = 0.148
EQ_soc PET 0.08 -0.12 0.27 p > .999
EQ_soc EYES 0.30 0.11 0.47 p = 0.066
PET EYES 0.11 -0.09 0.30 p > .999 
p3_t2 <- 
  ggplot(subscale_cors, aes(r)) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  geom_histogram(binwidth = 0.1, boundary = 0) +
  scale_x_continuous(limits = c(-1, 1), 
                     name = "Correlation", 
                     breaks = scales::breaks_pretty(n = 10)) + 
  scale_y_continuous(name = "Frequency", 
                     breaks = scales::breaks_pretty(n = 7)) + 
  theme_linedraw() +
  theme(panel.grid.minor.y = element_blank())

#p3_t2

# code repeated in order to use the rename, which screws up the table but is needed for the plot
p4_t2 <- dat |>
  select(SITES, 
         IRI_FS, IRI_EC, IRI_PT, IRI_PD, 
         EQ_cog, EQ_aff, EQ_soc, 
         PET_total,
         EYES_total) |>
  rename_with(~ str_remove(.x, "_total")) |>
  rename_with(~ str_replace(.x, "_", "\n")) %>%
  select(order(colnames(.))) |>
  correlation() %>%
  mutate(across(where(is.numeric), ~ round_half_up(., 2))) |>
  summary(redundant = TRUE) |>
  plot(show_data = "points",
       stars = FALSE) +
  theme_linedraw() +
  ggtitle("")

#p4_t2
p3_t2 + ggtitle("Correlation histogram") +
  p4_t2 + ggtitle("Correlation matrix") + 
  plot_layout(ncol = 1, heights = c(.4, 1))

ggsave("plots/cor_subscales_t2.png", 
       width = 7,
       height = 9,
       dpi = 600)

Both times - scales

p1_t1 + ggtitle("Correlation histogram (T1)") +
  p1_t2 + ggtitle("Correlation matrix (T2)") + 
  p2_t1 + ggtitle("Correlation histogram (T1)") +
  p2_t2 + ggtitle("Correlation matrix (T2)") + 
  plot_layout(ncol = 2, heights = c(.4, 1), guides = "collect")

ggsave("plots/cor_scales_both.png", 
       width = 9.35,
       height = 6,
       dpi = 600)

Both times - subscales

p3_t1 + ggtitle("Correlation histogram (T1)") +
  p3_t2 + ggtitle("Correlation matrix (T2)") + 
  p4_t1 + ggtitle("Correlation histogram (T1)") +
  p4_t2 + ggtitle("Correlation matrix (T2)") + 
  plot_layout(ncol = 2, heights = c(.4, 1), guides = "collect")

ggsave("plots/cor_subscales_both.png", 
       width = 13.2,
       height = 9,
       dpi = 600)

Session info

sessionInfo()
## R version 4.3.3 (2024-02-29)
## Platform: aarch64-apple-darwin20 (64-bit)
## Running under: macOS 15.4
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: Europe/Zurich
## tzcode source: internal
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] kableExtra_1.4.0     knitr_1.49           purrr_1.0.4         
##  [4] insight_0.19.10      see_0.8.3            correlation_0.8.4   
##  [7] mirt_1.40            lattice_0.22-6       scales_1.3.0        
## [10] ggrain_0.0.4         ggstance_0.3.7       patchwork_1.3.0.9000
## [13] ggplot2_3.5.2        janitor_2.2.1        lavaan_0.6-17       
## [16] psych_2.4.6.26       irr_0.84.1           lpSolve_5.6.20      
## [19] forcats_1.0.0        stringr_1.5.1        tidyr_1.3.1         
## [22] dplyr_1.1.4          readr_2.1.5         
## 
## loaded via a namespace (and not attached):
##  [1] mnormt_2.1.1         pbapply_1.7-2        polynom_1.4-1       
##  [4] gridExtra_2.3        permute_0.9-7        sandwich_3.1-1      
##  [7] rlang_1.1.6          magrittr_2.0.3       multcomp_1.4-28     
## [10] snakecase_0.11.1     compiler_4.3.3       mgcv_1.9-1          
## [13] systemfonts_1.0.6    vctrs_0.6.5          quadprog_1.5-8      
## [16] pkgconfig_2.0.3      crayon_1.5.3         fastmap_1.2.0       
## [19] labeling_0.4.3       pbivnorm_0.6.0       rmarkdown_2.29      
## [22] tzdb_0.5.0           ragg_1.3.0           bit_4.6.0           
## [25] xfun_0.49            cachem_1.1.0         jsonlite_1.8.9      
## [28] gghalves_0.1.4       Deriv_4.1.3          parallel_4.3.3      
## [31] cluster_2.1.6        R6_2.6.1             bslib_0.8.0         
## [34] stringi_1.8.7        lubridate_1.9.4      jquerylib_0.1.4     
## [37] estimability_1.5.1   Rcpp_1.0.14          zoo_1.8-14          
## [40] parameters_0.21.6    Matrix_1.6-5         splines_4.3.3       
## [43] timechange_0.3.0     tidyselect_1.2.1     rstudioapi_0.17.1   
## [46] yaml_2.3.10          vegan_2.6-4          codetools_0.2-20    
## [49] dcurver_0.9.2        tibble_3.2.1         withr_3.0.2         
## [52] bayestestR_0.13.2    coda_0.19-4.1        evaluate_1.0.1      
## [55] survival_3.5-8       ggpp_0.5.6           xml2_1.3.6          
## [58] pillar_1.10.2        generics_0.1.3       vroom_1.6.5         
## [61] hms_1.1.3            munsell_0.5.1        xtable_1.8-4        
## [64] glue_1.8.0           emmeans_1.10.2       tools_4.3.3         
## [67] mvtnorm_1.3-3        grid_4.3.3           datawizard_0.10.0   
## [70] colorspace_2.1-1     nlme_3.1-164         cli_3.6.4           
## [73] textshaping_0.3.7    viridisLite_0.4.2    svglite_2.1.3       
## [76] gtable_0.3.6         sass_0.4.9           digest_0.6.37       
## [79] GPArotation_2023.8-1 TH.data_1.1-3        farver_2.1.2        
## [82] htmltools_0.5.8.1    lifecycle_1.0.4      bit64_4.6.0-1       
## [85] MASS_7.3-60.0.1